/* ****************************************************************** **
**    OpenSees - Open System for Earthquake Engineering Simulation    **
**          Pacific Earthquake Engineering Research Center            **
**                                                                    **
**                                                                    **
** (C) Copyright 1999, The Regents of the University of California    **
** All Rights Reserved.                                               **
**                                                                    **
** Commercial use of this program without express permission of the   **
** University of California, Berkeley, is strictly prohibited.  See   **
** file 'COPYRIGHT'  in main directory for information on usage and   **
** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
**                                                                    **
** Developed by:                                                      **
**   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
**   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
**   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
**                                                                    **
** ****************************************************************** */
                                                                       
// $Revision: 1.6 $
// $Date: 2009-11-02 22:23:58 $
// $Source: /usr/local/cvs/OpenSees/EXAMPLES/Example1/main.cpp,v $


// File: ~/model/main.C
//
// Written: fmk 08/99
//
// Purpose: this file contains a C++ main procedure to perform the analysis
// of example1 (found in most documents). In the main() procedure:
//      1) each object of the domain, i.e. Nodes, Elements, Constraints,
//         and LoadPattern objects are created and then added to the Domain.
//      2) the components of the analysis object are constructed and then
//         the Analysis object is created.
//      3) the analysis is performed.
//      4) the results are printed - here the contents of Domain and end of
//         the analysis operation.

// standard C++ includes

#include <stdlib.h>

#include <OPS_Globals.h>
#include <StandardStream.h>

#include <ArrayOfTaggedObjects.h>

// includes for the domain classes
#include <Domain.h>
#include <Node.h>
#include <Truss.h>
#include <UniaxialMaterial.h>
#include <ElasticTimoshenkoBeam3d.h>
#include <ElasticBeam3d.h>
#include <ForceBeamColumn3d.h>
#include <ElasticMaterial.h>
#include <Steel02.h>
#include <SP_Constraint.h>
#include <MP_Constraint.h>
#include <LoadPattern.h>
#include <LinearSeries.h>
#include <NodalLoad.h>
#include <CrdTransf.h>
#include <PDeltaCrdTransf3d.h>
#include <CorotCrdTransf3d.h>
#include <Fiber.h>
#include <FiberSection3d.h>
#include <UniaxialFiber3d.h>
#include <LobattoBeamIntegration.h>
#include <UniformExcitation.h>
#include <GroundMotion.h>
#include <PathSeries.h>
#include <TimeSeriesIntegrator.h>
#include <TransientIntegrator.h>
#include <Newmark.h>
#include <VariableTimeStepDirectIntegrationAnalysis.h>

// includes for the analysis classes
#include <StaticAnalysis.h>
#include <AnalysisModel.h>
//#include <Linear.h>
#include <KrylovNewton.h>
#include <LineSearch.h>
#include <NewtonLineSearch.h>
#include <InitialInterpolatedLineSearch.h>
#include <TransformationConstraintHandler.h>
#include <DOF_Numberer.h>
#include <RCM.h>
#include <LoadControl.h>
#include <BandSPDLinSOE.h>
#include <BandSPDLinLapackSolver.h>
#include <SymSparseLinSOE.h>
#include <SymSparseLinSolver.h>
#include <CTestEnergyIncr.h>
#include <CTestNormDispIncr.h>
#include <ID.h>
#include <OPS_Stream.h>
#include <NodeRecorder.h>
#include <Recorder.h>
#include <math.h>

// init the global variabled defined in OPS_Globals.h
StandardStream sserr;
OPS_Stream *opserrPtr = &sserr;




// main routine
int main(int argc, char **argv)
{

    //  now create a domain and a modelbuilder
    //  and build the model
    // HSS Beam and column
    double BeamLength = 300;
    double BeamThick = 10;
    double BeamA;
    double BeamI;
    double BeamJ;
    double MG;
    double ElasticM = 205000.0;
    double PossionR = 0.3;

    BeamA = BeamLength * BeamLength - (BeamLength - 2 * BeamThick) * (BeamLength - 2 * BeamThick);
    BeamI = 1.0/12*(pow(BeamLength, 4) - pow((BeamLength - 2 * BeamThick), 4));
    BeamJ = pow(2, 0.5) * BeamI;
    MG = ElasticM/(2 * (1 + PossionR));

    Domain *theDomain = new Domain();
   
    // create the nodes using constructor:
    //          Node(tag, ndof, crd1, crd2)
    // and then add them to the domain
    Vector Trans1(3);
    Trans1(0) = 0;
    Trans1(1) = 0;
    Trans1(2) = 0;
    Vector Trans2(3);
    Trans2(0) = 0;
    Trans2(1) = 0;
    Trans2(2) = 0;
    Vector XZ(3);
    XZ(0) = 0;
    XZ(1) = 1;
    XZ(2) = 0;
    CrdTransf *Corot1 = new CorotCrdTransf3d(1, XZ, Trans1, Trans2);
    CrdTransf *PDelta1 = new PDeltaCrdTransf3d(2, XZ, Trans1, Trans2);

    // Create Nodes
    Node *node11 = new Node(1, 6, 0.0, 0.0, 0.0);
    Node *node12 = new Node(2, 6, 5000.0, 0.0, 0.0);
    Node *node13 = new Node(3, 6, 10000.0, 0.0, 0.0);
    Node *node21 = new Node(4, 6, 0.0, 0.0, 3000.0);
    Node *node22 = new Node(5, 6, 5000.0, 0.0, 3000.0);
    Node *node23 = new Node(6, 6, 10000.0, 0.0, 3000.0);
    Node *node31 = new Node(7, 6, 0.0, 0.0, 6000.0);
    Node *node32 = new Node(8, 6, 5000.0, 0.0, 6000.0);
    Node *node33 = new Node(9, 6, 10000.0, 0.0, 6000.0);

    Node *node14 = new Node(10, 6, 11000.0, 0.0, 0.0);
    Node *node15 = new Node(11, 6, 16000.0, 0.0, 0.0);
    Node *node24 = new Node(12, 6, 11000.0, 0.0, 3000.0);
    Node *node25 = new Node(13, 6, 16000.0, 0.0, 3000.0);
    Node *node34 = new Node(14, 6, 11000.0, 0.0, 6000.0);
    Node *node35 = new Node(15, 6, 16000.0, 0.0, 6000.0);

    // Set Nodal Mass
    Matrix nodemass(6, 6);
    nodemass(0, 0) = 1000;

    node11->setMass(nodemass);
    node12->setMass(nodemass);
    node13->setMass(nodemass);
    node14->setMass(nodemass);
    node15->setMass(nodemass);
    node21->setMass(nodemass);
    node22->setMass(nodemass);
    node23->setMass(nodemass);
    node24->setMass(nodemass);
    node25->setMass(nodemass);
    node31->setMass(nodemass);
    node32->setMass(nodemass);
    node33->setMass(nodemass);
    node34->setMass(nodemass);
    node35->setMass(nodemass);

    theDomain->addNode(node11);
    theDomain->addNode(node12);
    theDomain->addNode(node13);
    theDomain->addNode(node21);
    theDomain->addNode(node22);
    theDomain->addNode(node23);
    theDomain->addNode(node31);
    theDomain->addNode(node32);
    theDomain->addNode(node33);

    theDomain->addNode(node14);
    theDomain->addNode(node15);
    theDomain->addNode(node24);
    theDomain->addNode(node25);
    theDomain->addNode(node34);
    theDomain->addNode(node35);

    // create an elastic material using constriuctor:  
    //          ElasticMaterialModel(tag, E)

    UniaxialMaterial *Elastic1 = new ElasticMaterial(1, ElasticM);
    UniaxialMaterial *ElsPla1 = new Steel02(2, 345.0, ElasticM, 0.01, 19.0, 0.925, 0.15, 0, 1, 0, 1);
    UniaxialMaterial *RigidMat = new ElasticMaterial(3, ElasticM*100);

    ElasticBeam3d *Beam11 = new ElasticBeam3d(1, BeamA, ElasticM, MG, BeamJ, BeamI, BeamI, 4, 5, *PDelta1);
    ElasticBeam3d *Beam12 = new ElasticBeam3d(2, BeamA, ElasticM, MG, BeamJ, BeamI, BeamI, 5, 6, *PDelta1);
    ElasticBeam3d *Beam21 = new ElasticBeam3d(3, BeamA, ElasticM, MG, BeamJ, BeamI, BeamI, 7, 8, *PDelta1);
    ElasticBeam3d *Beam22 = new ElasticBeam3d(4, BeamA, ElasticM, MG, BeamJ, BeamI, BeamI, 8, 9, *PDelta1);
    ElasticBeam3d *Beam13 = new ElasticBeam3d(15, BeamA, ElasticM, MG, BeamJ, BeamI, BeamI, 12, 13, *PDelta1);
    ElasticBeam3d *Beam23 = new ElasticBeam3d(16, BeamA, ElasticM, MG, BeamJ, BeamI, BeamI, 14, 15, *PDelta1);

    theDomain->addElement(Beam11);
    theDomain->addElement(Beam12);
    theDomain->addElement(Beam21);
    theDomain->addElement(Beam22);
    theDomain->addElement(Beam13);
    theDomain->addElement(Beam23);

    int i = 1;
    int j = 1;
    double locx;
    double locy;
    double AFiber;
    Fiber *fibers[400];
    int tempno;
    for (i = 1; i < 11; i++)
    {
        for(j = 1; j < 11; j++)
        {
            locx = -150.0 + BeamLength / 10.0 * 0.5 + (i - 1) * BeamLength / 10.0;
            locy = 140.0 + BeamThick / 10.0 * 0.5 + (j - 1) * BeamThick / 10.0;
            AFiber = BeamThick*BeamLength/100.0;
            tempno = 10*(i - 1) + j - 1;
            fibers[tempno] = new UniaxialFiber3d(tempno + 1, *ElsPla1, AFiber, locx, locy);
        }
    }
    for (i = 1; i < 11; i++)
    {
        for(j = 1; j < 11; j++)
        {
            locx = -150.0 + BeamThick / 10.0 * 0.5 + (i - 1) * BeamThick / 10.0;
            locy = -140.0 + (BeamLength - 2*BeamThick) / 10.0 * 0.5 + (j - 1) * (BeamLength - 2*BeamThick) / 10.0;
            AFiber = BeamThick*(BeamLength - 2*BeamThick)/100.0;
            tempno = 10*(i - 1) + j - 1 + 100;
            fibers[tempno] = new UniaxialFiber3d(tempno + 1, *ElsPla1, AFiber, locx, locy);
        }
    }
    for (i = 1; i < 11; i++)
    {
        for(j = 1; j < 11; j++)
        {
            locx = -150.0 + BeamLength / 10.0 * 0.5 + (i - 1) * BeamLength / 10.0;
            locy = -150.0 + BeamThick / 10.0 * 0.5 + (j - 1) * BeamThick / 10.0;
            AFiber = BeamThick*BeamLength/100;
            tempno = 10*(i - 1) + j - 1 + 200;
            fibers[tempno] = new UniaxialFiber3d(tempno + 1, *ElsPla1, AFiber, locx, locy);
        }
    }
    for (i = 1; i < 11; i++)
    {
        for(j = 1; j < 11; j++)
        {
            locx = 140.0 + BeamThick / 10.0 * 0.5 + (i - 1) * BeamThick / 10.0;
            locy = -140.0 + (BeamLength - 2*BeamThick) / 10.0 * 0.5 + (j - 1) * (BeamLength - 2*BeamThick) / 10.0;
            AFiber = BeamThick*(BeamLength - 2*BeamThick)/100.0;
            tempno = 10*(i - 1) + j - 1 + 300;
            fibers[tempno] = new UniaxialFiber3d(tempno + 1, *ElsPla1, AFiber, locx, locy);
        }
    }

    SectionForceDeformation *ColumnSection[10];
    for (i = 0; i < 10; i++)
    ColumnSection[i] = new FiberSection3d(i + 1, 400, fibers);

    BeamIntegration *integrateon1 = new LobattoBeamIntegration();

    ForceBeamColumn3d *Column11 = new ForceBeamColumn3d(5, 1, 4, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column12 = new ForceBeamColumn3d(6, 4, 7, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column21 = new ForceBeamColumn3d(7, 2, 5, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column22 = new ForceBeamColumn3d(8, 5, 8, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column31 = new ForceBeamColumn3d(9, 3, 6, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column32 = new ForceBeamColumn3d(10, 6, 9, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column41 = new ForceBeamColumn3d(11, 10, 12, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column42 = new ForceBeamColumn3d(12, 12, 14, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column51 = new ForceBeamColumn3d(13, 11, 13, 10, ColumnSection, *integrateon1, *Corot1);
    ForceBeamColumn3d *Column52 = new ForceBeamColumn3d(14, 13, 15, 10, ColumnSection, *integrateon1, *Corot1);
    theDomain->addElement(Column11);
    theDomain->addElement(Column12);
    theDomain->addElement(Column21);
    theDomain->addElement(Column22);
    theDomain->addElement(Column31);
    theDomain->addElement(Column32);
    theDomain->addElement(Column41);
    theDomain->addElement(Column42);
    theDomain->addElement(Column51);
    theDomain->addElement(Column52);

    // create the single-point constraint objects using constructor:
    //          SP_Constraint(tag, nodeTag, dofID, value)
    // and then add them to the domain
    SP_Constraint *sp4 = new SP_Constraint(1, 0, 1);
    SP_Constraint *sp5 = new SP_Constraint(1, 1, 1);
    SP_Constraint *sp6 = new SP_Constraint(1, 2, 1);
    SP_Constraint *sp7 = new SP_Constraint(2, 0, 1);
    SP_Constraint *sp8 = new SP_Constraint(2, 1, 1);
    SP_Constraint *sp9 = new SP_Constraint(2, 2, 1);
    SP_Constraint *sp1 = new SP_Constraint(3, 0, 1);
    SP_Constraint *sp2 = new SP_Constraint(3, 1, 1);
    SP_Constraint *sp3 = new SP_Constraint(3, 2, 1);
    SP_Constraint *sp10 = new SP_Constraint(4, 0, 1);
    SP_Constraint *sp11 = new SP_Constraint(4, 1, 1);
    SP_Constraint *sp12 = new SP_Constraint(4, 2, 1);
    SP_Constraint *sp13 = new SP_Constraint(5, 0, 1);
    SP_Constraint *sp14 = new SP_Constraint(5, 1, 1);
    SP_Constraint *sp15 = new SP_Constraint(5, 2, 1);
    theDomain->addSP_Constraint(sp1);
    theDomain->addSP_Constraint(sp2);
    theDomain->addSP_Constraint(sp3);
    theDomain->addSP_Constraint(sp4);    
    theDomain->addSP_Constraint(sp5);    
    theDomain->addSP_Constraint(sp6);
    theDomain->addSP_Constraint(sp7);
    theDomain->addSP_Constraint(sp8);
    theDomain->addSP_Constraint(sp9);
    theDomain->addSP_Constraint(sp10);
    theDomain->addSP_Constraint(sp11);
    theDomain->addSP_Constraint(sp12);
    theDomain->addSP_Constraint(sp13);
    theDomain->addSP_Constraint(sp14);
    theDomain->addSP_Constraint(sp15);

    // Construct an ID for equalDOF
    ID vdofs(3);
    vdofs(0) = 0;
    vdofs(1) = 1;
    vdofs(2) = 2;

    Matrix mpconstr1(3, 3);
    for (i = 0; i < 3; i++)
    {
        mpconstr1(i, i) = 1;
    }

    MP_Constraint *mp1 = new MP_Constraint(6, 12, mpconstr1, vdofs, vdofs);
    MP_Constraint *mp2 = new MP_Constraint(9, 14, mpconstr1, vdofs, vdofs);
    theDomain->addMP_Constraint(mp1);
    theDomain->addMP_Constraint(mp2);

    // construct a linear time series object using constructor:
    //          LinearSeries()
   
    TimeSeries *theSeries = new LinearSeries();
   
    // construct a load pattren using constructor:
    //          LoadPattern(tag)
    // and then set it's TimeSeries and add it to the domain
   
    LoadPattern *theLoadPattern = new LoadPattern(1);
    theLoadPattern->setTimeSeries(theSeries);
    theDomain->addLoadPattern(theLoadPattern);
    opserr << "Finish\n";
    // construct a nodal load using constructor:
    //          NodalLoad(tag, nodeID, Vector &)
    // first construct a Vector of size 2 and set the values NOTE C INDEXING
    // then construct the load and add it to the domain

    Vector Node7Load(6);
    Node7Load(0) = 0.0;
    Node7Load(2) = -5000.0;
    Node7Load(1) = 0.0;
    Node7Load(3) = 0.0;
    Node7Load(4) = 0.0;
    Node7Load(5) = 0.0;
    Vector Node4Load(6);
    Node4Load(0) = 0.0;
    Node4Load(2) = -5000.0;
    Node4Load(1) = 0.0;
    Node4Load(3) = 0.0;
    Node4Load(4) = 0.0;
    Node4Load(5) = 0.0;
    Vector Node9Load(6);
    Node9Load(0) = 0.0;
    Node9Load(2) = -5000.0;
    Node9Load(1) = 0.0;
    Node9Load(3) = 0.0;
    Node9Load(4) = 0.0;
    Node9Load(5) = 0.0;
    Vector Node6Load(6);
    Node6Load(0) = 0.0;
    Node6Load(2) = -5000.0;
    Node6Load(1) = 0.0;
    Node6Load(3) = 0.0;
    Node6Load(4) = 0.0;
    Node6Load(5) = 0.0;
    Vector Node15Load(6);
    Node15Load(0) = 0.0;
    Node15Load(2) = -5000.0;
    Node15Load(1) = 0.0;
    Node15Load(3) = 0.0;
    Node15Load(4) = 0.0;
    Node15Load(5) = 0.0;
    NodalLoad *theLoad1 = new NodalLoad(1, 7, Node7Load);
    NodalLoad *theLoad2 = new NodalLoad(2, 4, Node4Load);
    NodalLoad *theLoad3 = new NodalLoad(3, 9, Node9Load);
    NodalLoad *theLoad4 = new NodalLoad(4, 6, Node6Load);
    NodalLoad *theLoad5 = new NodalLoad(5, 15, Node15Load);

    //addNodalLoad(NodelLoad *, int LoadPatternTag)
    theDomain->addNodalLoad(theLoad1, 1);
    theDomain->addNodalLoad(theLoad2, 1);
    theDomain->addNodalLoad(theLoad3, 1);
    theDomain->addNodalLoad(theLoad4, 1);
    theDomain->addNodalLoad(theLoad5, 1);

    AnalysisModel     *theModel = new AnalysisModel();
    StaticIntegrator  *theIntegrator = new LoadControl(0.05, 10, 0.05, 0.05);
    ConstraintHandler *theHandler = new TransformationConstraintHandler();
    RCM               *theRCM = new RCM();

    DOF_Numberer      *theNumberer = new DOF_Numberer(*theRCM);
    //BandSPDLinSolver  *theSolver = new BandSPDLinLapackSolver();
    //LinearSOE         *theSOE = new BandSPDLinSOE(*theSolver);
    SymSparseLinSolver  *theSolver = new SymSparseLinSolver();
    SymSparseLinSOE *theSOE = new SymSparseLinSOE(*theSolver, 1);
    ConvergenceTest   *theTest = new CTestEnergyIncr(0.000001, 200, 1, 2, 1);

    LineSearch *theSearch = new InitialInterpolatedLineSearch();
    EquiSolnAlgo      *theSolnAlgo = new NewtonLineSearch(*theTest, theSearch);
    StaticAnalysis    theAnalysis(*theDomain, *theHandler, *theNumberer, *theModel, *theSolnAlgo, *theSOE, *theIntegrator, theTest);

    // perform the analysis & print out the results for the domain
    int numSteps = 50;
    int nstep = 0;
    for (nstep = 0; nstep < numSteps; nstep++)
    {
        theAnalysis.analyze(1);
    }
    StaticIntegrator  *theIntegrator2 = new LoadControl(0.01, 10, 0.01, 0.01);
    StaticAnalysis theAnalysis2(*theDomain, *theHandler, *theNumberer, *theModel, *theSolnAlgo, *theSOE, *theIntegrator2, theTest);
    opserr << "Changing.................\n\n";
    for (nstep = 0; nstep < 2; nstep++)
    {
        theAnalysis2.analyze(1);
    }

    // Start Triansient Analysis
    // PathSeries(int tag, const char *fileName, double pathTimeIncr = 1.0, double cfactor = 1.0, bool useLast = false, bool prependZero = false, double startTime = 0.0);
    TimeSeries *theSeismicWave = new PathSeries(2, "/home/chenting/Chenting/TestExp/OpenSees/EXAMPLES/Example2/NGA01.txt", 0.01, 9.8, true, true, 0);
    theDomain->setCurrentTime(0.0);
    LoadPattern *theSeismicPattern = new LoadPattern(2);
    theSeismicPattern->setTimeSeries(theSeismicWave);
    theDomain->addLoadPattern(theSeismicPattern);

    TimeSeries *EmptyVelSeries = 0;
    TimeSeries *EmptyDispSeries = 0;
    TimeSeriesIntegrator *EmptyIntegrator = 0;

    //GroundMotion *GM = new GroundMotion(dispSeries,velSeries,accelSeries,seriesIntegrator,dtInt,fact);
    //It's useless to modify dtInt and fact
    GroundMotion *GM = new GroundMotion(EmptyDispSeries, EmptyVelSeries, theSeismicWave, EmptyIntegrator, 0.01, 1.0);

    //UniformExcitation(GroundMotion &theMotion, int dof, int tag, double vel0 = 0.0, double fact = 1.0);
    UniformExcitation *theExcitation = new UniformExcitation(*GM, 1, 1);

    // Define ImposedMotion

    TransientIntegrator *theIntegrator3 = new Newmark(0.5, 0.25);

    VariableTimeStepDirectIntegrationAnalysis theAnalysis3(*theDomain, *theHandler, *theNumberer, *theModel, *theSolnAlgo, *theSOE, *theIntegrator3, theTest);
    for (nstep = 0; nstep < 100; nstep++)
    {
        theAnalysis3.analyze(1, 0.01, 0.0001, 0.05, 20);
        opserr << *theDomain;
    }

    //theAnalysis3.analyze(2000, 0.01, 0.0001, 0.05, 20);

    exit(0);
}